FILE_LOCATION   =	"/ifs/scratch/c2b2/ip_lab/ddd2124"

setwd(FILE_LOCATION)

HSIB_FILE<-read.table("halfsiblings.lst",header=F,stringsAsFactors=F)

poplist<-read.table("PopIndex.txt",header=F, stringsAsFactors=F)
poplist<-poplist[,1]



#enact these once you have the lengths and numbers files
Lengths<-read.table("Lengths.txt",header=F)
Numbers<-read.table("Numbers.txt",header=F)



hsibl<-integer(0)
hsibn<-integer(0)

for (i in seq(length(HSIB_FILE[,1])))
{
	if ((A<-which(poplist==HSIB_FILE[i,1]))<(B<-which(poplist==HSIB_FILE[i,2])))
	{hsibl<-c(hsibl,Lengths[A,B])
	 hsibn<-c(hsibn,Numbers[A,B])}
	else 	{hsibl<-c(hsibl,Lengths[B,A])
		 hsibn<-c(hsibn,Numbers[B,A])}
}

hsiblcdf<-ecdf(hsibl)
hsibncdf<-ecdf(hsibn)


g<-length(Lengths[,1])

#####################################################
##################HSIB MATRICES######################

HSIBNMAT <- Lengths
for (i in seq(g-1))
{
	HSIBNMAT[i,] <- hsibncdf(Numbers[i,])
	HSIBNMAT[i,] <- apply(rbind(HSIBNMAT[i,],1-HSIBNMAT[i,]), 2, min)
}

HSIBLMAT <- PCNMAT
for (i in seq(g-1))
{
	HSIBLMAT[i,] <- hsiblcdf(Lengths[i,])
	HSIBLMAT[i,] <- apply(rbind(HSIBLMAT[i,],1-HSIBLMAT[i,]), 2, min)
}

HSIBNMAT2 <- HSIBNMAT + t(HSIBNMAT)
write.table(HSIBNMAT2,"HSIBNMAT.txt",row.names=F,col.names=F,quote=F)
rm(HSIBNMAT2, HSIBNMAT)
gc()
HSIBLMAT2 <- HSIBLMAT + t(HSIBLMAT)
write.table(HSIBLMAT2,"HSIBLMAT.txt",row.names=F,col.names=F,quote=F)
rm(HSIBLMAT2, HSIBLMAT)
gc()
###############################################################################
# add matrices for avg lengths

AverageLengths <- Numbers / (Lengths + .001)
hsiba<-integer(0)

g<-length(AverageLengths[,1])

for (i in seq(length(HSIB_FILE[,1])))
{
	if ((A<-which(poplist==HSIB_FILE[i,1]))<(B<-which(poplist==HSIB_FILE[i,2])))
	{hsiba<-c(hsiba,AverageLengths[A,B])}
	else 	{hsiba<-c(hsiba,AverageLengths[B,A])}
}

hsibacdf<-ecdf(hsiba)

HSIBAMAT<- AverageLengths
for (i in seq(g-1))
{
	HSIBAMAT[i,] <- hsibacdf(AverageLengths[i,])
	HSIBAMAT[i,] <- apply(rbind(HSIBAMAT[i,],1-HSIBAMAT[i,]), 2, min)
}

HSIBAMAT<- HSIBAMAT+ t(HSIBAMAT)
write.table(HSIBAMAT,"HSIBAMAT.txt",row.names=F,col.names=F,quote=F)
rm(HSIBAMAT, hsiba, hsibacdf)
gc()

###################

quit()